﻿/*
 * StandardMolodenskyTransform
 * 
 * initial implementation : 20.11.2008
 *
 * Copyright (c) Adam Sida, asida@centrum.cz.
 * All rights reserved.
 * 
 * for license see Petra.Gis.License.txt
 */
using System;
using System.Collections.Generic;
using System.Runtime.Serialization;
using Petra.Gis.Geometry.Support;


namespace Petra.Gis.ReferenceSystemServices.Transformations
{

  /// <summary>
  /// StandardMolodensky is geographic space transform defined
  /// as simple 3D shift between source and target spheroid.
  /// </summary>
  [DataContract]
  public class StandardMolodenskyTransform : DatumTransform
  {

    private double from_a;
    private double from_f;
    private double from_esq;

    private double da;
    private double df;

    // transform variables
    /// <summary>
    /// Molodensky shift vector in meters.
    /// </summary>
    private Vector pos;

    //semi major / semi minor
    private double adb;
    double slat, clat, slon, clon, ssqlat;
    double dlat, dlon, dh;
    double rn, rm;
    //double from_Z;

    Coordinate out3D;




    #region constructors


    /// <summary>
    /// xaml constr.
    /// </summary>
    public StandardMolodenskyTransform() 
      :base()
    {
      Metadata = new Metadata("StandardMolodenskyTransform", "none", -1, "", "", "provides transformation between 2 Geographics Coordinate systems with different spheriods");
    }


		public StandardMolodenskyTransform(GeographicCoordinateSystem sourceCS, 
                                       GeographicCoordinateSystem targetCS,
                                       DistanceUnit distanceParamUnit,
                                       double dx, double dy, double dz, bool createInverse)
			:this(sourceCS, targetCS, distanceParamUnit, createInverse,
        new XParam("TranslationX", dx),
        new XParam("TranslationY", dy),
        new XParam("TranslationZ", dz)
      )
    {  }



    

    public StandardMolodenskyTransform(GeographicCoordinateSystem sourceCS,
                                       GeographicCoordinateSystem targetCS,
                                       DistanceUnit distanceParamUnit, bool createInverse, params XParam[] parameters)
      : this(sourceCS, targetCS, distanceParamUnit, createInverse, new TransformParameters(parameters) )
    {  }

    
    
    
    
    
    public StandardMolodenskyTransform(GeographicCoordinateSystem sourceCS,
                                       GeographicCoordinateSystem targetCS,
                                       DistanceUnit distanceParamUnit, bool createInverse, TransformParameters parameters)
      : base(sourceCS, targetCS, distanceParamUnit, null, createInverse, parameters)
    {
      if (parameters.Count != 3)
        throw new ArgumentException("StandardMolodenskyTransform requires 3 parematers");

      Metadata = new Metadata("StandardMolodenskyTransform", "none", -1, "", "", "provides transformation between 2 Geographics Coordinate systems with different spheriods");
    }


    #endregion constructors


    //protected override void transform(ICoordinate srcCoord)
    //{
    //  transform(srcCoord, Dimension.dim3D);
    //}




    protected override void transformN(ICoordinate from)
    {

      // convert source units to radians
      //out3D.lon = from.lon * srcUnit.RadiansPerUnit;
      //out3D.lat = from.lat * srcUnit.RadiansPerUnit;

      //if (dataDim == Dimension.dim2D)
      //  out3D.Z = 0;
      //else
      out3D.Z = from.Z;

      slat = Math.Sin(out3D.lat); //latitude
      clat = Math.Cos(out3D.lat);
      slon = Math.Sin(out3D.lon); //longitude
      clon = Math.Cos(out3D.lon);

      ssqlat = slat * slat;

      rn = from_a / Math.Sqrt(1.0 - from_esq * ssqlat);
      rm = from_a * (1.0 - from_esq) / Math.Pow((1.0 - from_esq * ssqlat), 1.5);

      dlat = (((((-pos.I * slat * clon - pos.J * slat * slon) + pos.K * clat)
                  + (da * ((rn * from_esq * slat * clat) / from_a)))
              + (df * (rm * adb + rn / adb) * slat * clat)))
          / (rm + out3D.Z);

      dlon = (-pos.I * slon + pos.J * clon) / ((rn + out3D.Z) * clat);

      dh = (pos.I * clat * clon) +
           (pos.J * clat * slon) +
           (pos.K * slat) -
           (da * (from_a / rn)) +
           ((df * rn * ssqlat) / adb);

      out3D.lon += dlon;
      out3D.lat += dlat;
      out3D.Z += dh;

      // convert to target angular unit
      //from.lon /= trgtUnit.RadiansPerUnit;
      //from.lat /= trgtUnit.RadiansPerUnit;

      //if (dataDim == Dimension.dim3D)
      from.Z = from.Z;

    }





    protected override bool parametersValid()
    {
      return (Parameters.Count == 3);
    }




    /// <summary>
    /// must override here, spheric conformal doesnt need parameter angular units (only distance)
    /// </summary>
    /// <returns></returns>
    protected override bool graphCompleted()
    {
      return !(this.DistanceParameterUnit == null || !parametersValid());
    }





    protected override void updateTransform()
    {

      // prepare spheroid parameters
      GeographicCoordinateSystem cs1, cs2;
      if (!IsInverse)
      {
        cs1 = (GeographicCoordinateSystem)Source;
        cs2 = (GeographicCoordinateSystem)Target;
      }
      else
      {
        cs1 = (GeographicCoordinateSystem)Target;
        cs2 = (GeographicCoordinateSystem)Source;
      }

      from_a = cs1.Datum.Ellipsoid.SemiMajorAxis * cs1.Datum.Ellipsoid.AxisUnit.MetersPerUnit;
      from_f = cs1.Datum.Ellipsoid.InverseFlattening;
      from_esq = Math.Pow(cs1.Datum.Ellipsoid.Eccentricity, 2.0);

      double to_a = cs2.Datum.Ellipsoid.SemiMajorAxis * cs2.Datum.Ellipsoid.AxisUnit.MetersPerUnit;

      da = to_a - from_a;
      df = cs2.Datum.Ellipsoid.InverseFlattening - cs1.Datum.Ellipsoid.InverseFlattening;

      adb = 1.0 / (1.0 - from_f);  // "a divided by b"

      pos = new Vector();

      pos.I = Parameters["TranslationX"];
      pos.J = Parameters["TranslationY"];
      pos.K = Parameters["TranslationZ"];

      // normalize intput parameters units
      pos.I *= base.DistanceParameterUnit.MetersPerUnit;
      pos.J *= base.DistanceParameterUnit.MetersPerUnit;
      pos.K *= base.DistanceParameterUnit.MetersPerUnit;

      //for inverse must 
      //change sign for pos, rot and m sign
      if (IsInverse)
      {
        pos.I *= -1;
        pos.J *= -1;
        pos.K *= -1;
      }

      out3D = new Coordinate();
    }




    public override SRSTransform Inverse()
    {
      return new StandardMolodenskyTransform(
                      (GeographicCoordinateSystem)Source,
                      (GeographicCoordinateSystem)Target,
                      base.DistanceParameterUnit,
                      !IsInverse,
                      Parameters
                      );
    }





    public override Dimension TransformDimension
    {
      get { return Dimension.dim3D; }
    }


  }
}
